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ABSTRACT 



We present a new approach to calculate the Wiener filter solution of general data sets. It is trivial to implement, flexible, 
numerically absolutely stable, and guaranteed to converge. Most importantly, it does not require an ingenious choice 
of preconditioner to work well. The method is capable of taking into account inhomogeneous noise distributions and 
arbitrary mask geometries. It iterative ly builds up the signal reconstruction by means of a messenger field, introduced 
to mediate between the different preferred bases in which signal and noise properties can be specified most conveniently. 
Using cosmic microwave background (CMB) radiation data as a showcase, we demonstrate the capabilities of our scheme 
by computing Wiener filtered WMAP7 temperature and polarization maps at full resolution for the first time. We show 
how the algorithm can be modified to synthesize fluctuation maps, which, combined with the Wiener filter solution, 
result in unbiased constrained signal realizations, consistent with the observations. The algorithm performs well even 
on simulated CMB maps with Planck resolution and dynamic range. 
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1. Introduction 

Signal reconstruction out of noisy data is arguably one of 
the most fundamental problems in astronomy (and science 
in general), and has been subject of extensive research for 
centuries. The method of least squares, for example, was 
originally introduced to predict the orbital parameters of 
the dwarf p lanet 1 Cere s out of a series of astrometric mea- 
surements (|Gauss|[l809l ). Since then, applications in signal 
processing have become legion. 

Incorporating statistical information about the proper- 
ties of noise and signal laid the foundation for a large class 
of powerful reconstruction techniques, including particu- 
larly popular methods such as the Wiener filter ([Wiener! 
194 9]), appro aches that rely on entropy maxi mization 
([Javnedfl957l h or the Kalman filter (|Kalmanlll960h . The is- 
sue has remained a topic of ongoing research. For Gaussian 
fields, only recently efficient algorithms have been de- 
veloped for case s without prior knowledge of the signal 
properties (e.g. , Wandelt et all 12004 Uasche et all 120101: 
lEnfflin fc Frommertl l20lil T or, where neith er signal nor 
noise covariances are known a priori ([Oppermann et all 

unnt). 

In this letter, we focus on the evaluation of the gen- 
eralized Wiener filter. If the observed data d are a linear 
combination of signal s and noise n, 

d = s + n , (1) 

the Wiener filter is defined as the solution of the equation 

(S- 1 +N- 1 )s WF =N- 1 d. (2) 

Then, by construction, swf is the maximum a posteriori 
probability solution if signal and noise are both Gaussian 



random fields with known covariances S and N, respec- 
tively. 

The Wiener filter or variants thereof are fre- 
quently encountered in the post-processing of observa- 
tional data. In the analysis of cosmic microwave back- 
ground (CMB) radiation experiments, for exam ple, dur- 
ing optimal power spe ctrum estimati on (e.g., [ Tegmark 
1997bt iBond et al.lll998l lOh et al.l ll 999t lElsner fc Wandeltl 
201 2a|), likelihood analysis (e.g., [H inshaw jet al.l 120071 : 
iDunklev etafl2009t [ Eisner fc Wandeltll2012bD . mapmaking 
(e.g. jTegmark et al.lll997HTegmarklll997a|), and lensing re - 
constructions (e.g., iHirata et al.l 12004 ISmith et al.l 120071 ) , 
etc. 

In practice, the numerical computation of Eq. (j2]) can be 
challenging for the large data sets collected by modern ex- 
periments. The underlying difficulty is that signal and noise 
covariances grow as the square of the size of the data set and 
become too large to be stored and processed as dense matri- 
ces. Progress is possible when these matrices are structured, 
e.g. when there are sets of bases where S and N become 
sparse. However, the bases may be incompatible, for exam- 
ple, the signal covariance may be sparse in Fourier space, 
whereas the noise covariance may be sparse in pixel space. 
As a result, one can either solve the Wiener filter equation 
only approximately, e.g. by assuming a hom ogeneous and 
isotropic noise distribution fas done in, e.g., IHirata et al.l 
120041 : iKomatsu et~aIll2QQ5l : iMangilli fc Verdell2009l ). or has 
to rely on complex and costly numerical schemes (e.g., in- 
volving c onjugate gradient solvers with multigrid precon- 
ditioner s, |SmitF^t_^i] [2003) . Finding fast preconditioners 
that work is a highly non-trivial art especially since the ma- 
trices of interest are often extremely ill-conditioned, with 
typical condition numbers of the order (D(l0 7 ). 
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The Wiener filter is the optimal linear filter only for 
strictly Gaussian noise and signal. However, a maximum 
entropy argument shows that the Gaussian prior is the 
least informative prior, once mean and covariance of the 
random fields are specified. From a Bayesian point of view, 
a Gaussian prior can therefore be a well-justified approx- 
imation even for non- Gaussian inference. This argument 
provides a theoretic al underpinning for the work of, e.g, 
Bou chet et all ([1999), where the authors use Wiener filter- 
ing as a means to clean CMB maps from (non-Gaussian) 
foregrounds. Applied in the limit of weak non-Gaussianity, 
the Wiener filter has also proven an indispensable tool 
for th e search for primor dial non-Gaussianity in CMB 
data (|Komatsu et al.l 12005), where inference from the full 
non- Gaussian posterior is computationally very demanding 
(lElsner fc Wandeltll2010l) . 

In this letter, we propose a high precision, iterative al- 
gorithm for the solution of the exact Wiener filter equation. 
It is conceptually very simple and therefore easy to imple- 
ment correctly, yet still numerically efficient enough to be 
of general interest. The method is fully universal as long as 
there exist easily accessible orthogonal basis sets which ren- 
der S and N sparse. We illustrate our approach by applying 
it to CMB signal reconstruction. 

The letter is organized as follows. In Sect. [2j we intro- 
duce our new algorithm for the efficient calculation of the 
Wiener filter. We then present an explicit example and ana- 
lyze CMB data of the WMAP experiment (Sect. [3} showing 
how to generalize Eq. ([T]) to include an instrumental trans- 
fer function. Finally, we summarize our findings in Sect. SJ 
Throughout the paper, we assume WMAP7+ BAO+H0 cos- 
mological parameters (Ko matsu et al.|[201lf ). 

2. Method 

In the following, we present the basic equations of the al- 
gorithm and address how to solve them numerically before 
discussing a concrete example. 

2.1. Algorithm for Wiener filtering 

For an efficient evaluation of Eq. (j2j), we introduce an ad- 
ditional stochastic auxiliary field t with covariance T. The 
statistical properties of this helper variable are deliberately 
chosen very simple: T is proportional to the identity ma- 
trix. The identity matrix is invariant under any orthogonal 
basis change. While it may not be possible to directly apply 
expressions like (S + N) _1 to a data vector, we can always 
do so with combinations as (S + T) _1 and (N + T) _1 re- 
gardless of which basis renders S or N sparse. 

In order to take advantage of this possibility, we form 
a set of equations to be solved iteratively for the signal 
reconstruction s and the auxiliary field t, 

(INT 1 + (AT)" 1 ) t = INT 1 d + (AT)" 1 s (3) 
(S- 1 + (AT)- 1 ) s = (\T)-H, (4) 

where we defined N = N — T and introduced one additional 
scalar parameter A, whose usefulness we will discuss in the 
context of convergence acceleration in Sect. 12.21 We choose 
the amplitude of the covariance matrix of the auxiliary field 
according to T = rt with r = min(diag(N)). Consequently, 
the field t gains a physical interpretation: it corresponds to 



a signal reconstruction given the fraction of the noise that 
is distributed homogeneously and isotropically in the data. 
It can be easily shown that the above system reduces to the 
Wiener filter equation in the limit A = 1. 

The algorithm can be summarized as follows. We initial- 
ize the vectors s and t with zeros. First, we solve Eq. (|3]) 
for the auxiliary field t in the basis defined by the noise co- 
variance matrix. Next, we change the basis to, e.g., Fourier 
space, where S can be represented easily. Then, we solve 
for s using Eq. (jl]) given the latest realization of t, and fi- 
nally transform the result back to the original basis. After 
a sufficient number of iterations, the signal reconstruction 
converges to the Wiener filter solution, s — ^ swf as we show 
now. 

2.2. Stability, convergence speed and diagnostics 

The method is unconditionally stable. If the residual is 
at one iteration, it will be 

e i+1 = [S(S + AT)" 1 ] [N(N + AT)" 1 ] a (5) 

at the next. If N ex 1, the system is solved in a single 
iteration. The terms in square brackets have spectral radii 
pi 5 2 strictly less than unity. So |e^ + i| < P1P2 |^| for all 
z, and convergence is exponential. It is fast for low noise 
pixels and signal modes with low prior variance, and slowest 
for modes with high prior variance in high noise pixels. 
Temporarily setting A > 1 accelerates convergence in this 
problematic regime. 

For a given A, the iterations minimize 

X 2 (x, A) = x t S- 1 x+(^-x) t [N + (A - 1)T] _1 (d-x) (6) 

to give s(X) = S(S + N + (A - 1)T) _1 . Fortunately, it is 
precisely in the regime where convergence is slowest (S, N 
large) that s(X > 1) approximates best. 

The goal is therefore to find a "cooling schedule" which 
starts at high A and reduces it to A = 1 such that the 
final solution satisfies a defined accuracy criterion. Note 
that further iterations at lower A never degrade but always 
improve parts of the solution that were found at higher A. 

At x = s(A), we have Xmin( A ) = ^ f ( s + N + (A - 
^T)- 1 ^ = dtS-yA) < XminW for A > 1. We now show 
that we can find a "cooling schedule" for A which respects 
the error bound E = x 2 (A,x) < b = XmmiX) throughout. 

The following algorithm has the desired property: ini- 
tialize x = and choose the initial A such that x 2 (A,0) = 
d f (N + (A - l)T)- 1 (i = b. Then iterate x once at that A, 
which reduces E to x 2 (A,x). Decrease A such that E = b. 
Repeat until A = 1 at which point x = s(l). 

So far we have not specified how to calculate the opti- 
mal b = Xmin(l)> wm ch only becomes easy to compute after 
the Wiener Filter has been found. Choosing b > Xmin(l) 
leads to a needlessly inaccurate map, containing numeri- 
cal noise artifacts that could result in erroneous statistical 
conclusions; choosing b < XminW? gi yes results that are not 
wrong, just suboptimal. 

Fortunately, it is possible to approximate b during the 
solution, as long as the initial b < XminW- After converging 
at A > 1 when b < XminW? we can increase b to make 
it sufficiently close to XminW m one ste P w ^h minimal 
additional computation. At that point, we have x = s(A) 
and can calculate b' = Xmin(^) + — l)rx^S~ 2 x < XminW? 
with the inequality being quadratically small. 
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If d is consistent with the assumed statistical model, 
d — s + n, we can find how close b should be to XminW- 
For an ensemble of data sets, XminW * s X 2_ distributed with 
a mean given by the number of degrees of freedom, na.o.fo 
and a variance of <r^ 2 = 2na.o.f.- Hence, no outcome of 
statistical inference based on the resulting Wiener Filter 
depends on whether it has an excess x 2 (l) °f much less 
than <7 % 2. 

Even if we had chosen b far below the expectation, 
b = Ttd.o.f — rncr x 2, with m > 5, excluding a violation of 
b < Xmin(l) f° r an " Practical purposes, the updated b' satisi- 
fies Xmin(l) — b'~ a x 2/ v /n c i.o.f.- It is therefore statistically 
irrelevant as long as nd.o.f. ^ 1 which is true for all cases 
where an iterative scheme is likely to be useful. 

This algorithm guarantees satisfying a given error 
bound, but it does not guarantee the fastest possible time 
to solution. For applications to the CMB, where the signal 
power spectrum tends to decrease with increasing scale, it is 
much cheaper to iterate at high A than at low A. In Sect.[3j 
we describe a heuristic choice for the "cooling schedule" 
that exploits the assumed shape of the signal power spec- 
trum and leads to significantly faster execution times. 

2.3. Constrained realizations 

The Wiener filter is known to remove noise very aggres- 
sively. As a result, the covariance of the signal reconstruc- 
tion is reduced compared to the true signal, 

(s W f4f) = S(S + N)- 1 S. (7) 

To obtain signal realizations with correct covariance prop- 
erties that are consistent with the observed data, we have 
to add a fluctuation v ector / to the W iener filter so- 
lution with covariance ([Bertschingerl Il987l see also, e.g., 
iHoffman fc Ribaklll99l[ ) 

(// t ) = (S- 1 +N- 1 )- 1 . (8) 

Applying only minor modifications, the algorithm as intro- 
duced above can be adapted to generate suitable fluctua- 
tion maps. The only difference lies in the source term on 
the right-hand side of Eq. (j3j) , which we replace by 

1ST 1 d -> N _1 N (l\T 1/2 gi + S~ 1/2 2 ) , (9) 

where gi^ are independent univariate Gaussian random 
variables. For those entries of the vector S -1//2 #2 that fall 
within masked regions where N _1 is zero and formally no 
inverse N can be defined, the corresponding blocks of the 
matrix product N _1 N have to be replaced with the identity 
matrix. 

After the algorithm has been evaluated as explained 
in Sect. 12. 1[ the solution is a random realization of a 
fluctuation map / with the correct covariance properties. 
Constrained realizations are then given by sqr = swf + /• 

3. WMAP polarization analysis 

We now discuss the application of our technique to 
Wilkinson microwa ve anisotropy probe ( WMAP) data 
([Jarosik et al.ll20ilh . We used the seven-year observations 



of the CMB radiation in the V-band at 61 GHz as a ba- 
sis for the reconstruction^], at resol ution parameter n^ de = 
512 in the Healpix representation ([Gorski et al.ll2QQ5[ ) and 
f max = 1024. Since CMB anisotropics are isotropic and 
Gaussian, the signal covariance is diagonal in spherical 
harmonic space. The noise covariance, on the other hand, 
can be approximated by a diagonal matrix in pixel space 
(jffinshaw et al.ll2003D . 

However, as we are interested in the Wiener filtering 
of temperature and polarization data, S and N are in 
fact block diagonal, with a si gnificantly increase d condi- 
tion number of the matrices ([Larson et al.l [20071 ). For all 
multipole moments I, we write, 

/^TT /nrTE\ 

St = b " \c™ cfV ' (10) 

where the Cf- X are the power spectra for temperature (XX 
= TT), polarization (XX = EE), and their cross-spectrum 
(XX = TE), and the bg are the beam function. As we did 
not aim at a reconstruction of a possible contribution from 
a curl component to the polarization signal (XX = BB), 
we excluded it from the analysis. The blocks of the noise 
covariance matrix are, for every pixel, 

/(II) \ 
Ni= (QQ) (QU) , (11) 
V (QU) (UU)/ 

where we used the Stokes parameters I, Q, and U. For 
WMA P data, the cro ss-correlations (IQ) and (IU) are 
small ([Page et al.ll2007D . 

We adopted the official WMAP extended temperature 
and polarization masks, restricting the sky fraction avail- 
able to 70.6 % and 73.3 %, respectively. We removed pos- 
sible residual monopole and dipole contributions after the 
masks had been applied. Avoiding numerical errors asso- 
ciated with spherical harmonic transforms on the cut sky 
can be done by an explicit low rank update of the inverse 
noise covariance matrix N _1 in Eq. (J3j) to assign infinite 
variance to a set of template maps, rep resenting the (non - 
cosmological) modes to be subtracted ([Smith et al.ll2009l ). 
However, we found it numerically more stable to split the 
calculation into separate steps; we therefore multiplied the 
input independently by 

d= lim (t + aVV^y 1 d raw , (12) 

where V is a 4 x n p i x matrix, storing the four masked tem- 
plates to be marginalized over. We used the Woodbury ma- 
trix identity to evaluate the equation. 

To initialize the algorithm, we chose A according to 
the ratio T/Sj^ at i^er — 20. This starting point is well 
suited to let all large-scale modes from t w 30 down to the 
quadrupole converge. After the x 2 nas reached a plateau, we 
decreased it to the ratio at if™ = 2^ e d r . We restricted the 
maximal increment in £i ter to 500 and assured that A > 1 
is always fulfilled. We achieved robust results after termi- 
nating the algorithm once A\ 2 < 10 -4 <7 % 2 at A = 1 has 
been reached, in this case at a final x 2 / n d.o.f. = 0.996. The 
success of the algorithm does not sensitively depend on the 
scheme for lowering A to unity. 

1 Available at http : //lambda . gsf c . nasa . gov 
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The gradual decrease of A, coupled to a well specified 
multipole moment ^i te r, allows for a significant speedup of 
the algorithm. Since the signal is build-up from large to 
small scales, fluctuations at multipoles i > l- lteT will be 
suppressed. As a consequence, we can restrict the max- 
imum multipole moment of the computationally expen- 
sive spherical harmonic transforms to ~ iiter + 100, 
considerably smaller than ^ max for most of the iterations. 
Owing to a time complexity of (9(n sMe ^max) , the over- 
all computing time for the spherical harmonic transforms 
can be substantially reduced. For the final iterations in- 
cluding the full range of the algorithm is an excellent 
candidate for further acceleration by means of the mas- 
sive parallelism offered by modern graphics processing units 
(jElsner fc Wandeltll201lh . 

The result of the WMAP analysis is shown in Figs. [1] 
and where we plot patches of the Wiener filtered maps 
for the Stokes parameters I, Q, and U, as well as an exam- 
ple of a temperature fluctuation map. The non- vanishing 
cross-spectrum Cj E allows for an efficient polarization re- 
construction in masked regions, as long as temperature data 
are available. A distinguishing feature of Wiener filtering is 
the extrapolation of signal into the mask; this works es- 
pecially well in the small areas where single point sources 
have been cut out. 

4. Discussion and conclusions 

We presented an efficient technique to calculate the Wiener 
filter solution of general data sets. Introducing a stochastic 
auxiliary field with a simple covariance matrix proportional 
to the identity matrix, we indicated how an iterative solver 
for the Wiener filter equation can be constructed. The re- 
sulting algorithm is not only numerically robust, but also 
flexible. Using CMB temperature and polarization data ob- 
tained by the WMAP mission as an example, we explicitly 
demonstrated the usefulness of the method. 

We estimated the numerical efficiency of the algorithm 
for typical CMB temperature maps as obtained with the 
high frequency instrument aboa rd the Planck satellite 
([Planck Collaboration et all l201lf ) . At Healpix resolution 
^side = 2048 and ^ max = 3000, calculating the Wiener filter 
solution takes about 22 CPUhrs on a single computing node 
comprising two quad-core Intel Xeon 2.7 GHz processors. 
We conclude, that the method described here is not only 
applicable to a wide range of problems, but also effective 
enough to be of general interest. We are currently exploring 
several further variants of the general messenger approach 
as well as generalizations of this method to the case of mul- 
tiple data sets measuring the same underlying signal with 
different transfer functions or more general noise covariance 
matrices. 
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Fig. 1. Temperature reconstruction. Using the WMAP mask (left panel) for the analysis, we plot the V-band temperature 
map after Wiener filtering (middle panel). It can be supplemented to a constrained realization by superposing a fluctuation 
term (right panel). The patches are 25° on the side, centered at galactic coordinates (I, b) = (0°, 35°). 




Fig. 2. Polarization reconstruction. Using the WMAP polarization mask (left panel) for the analysis, we obtain Wiener 
filtered maps for the Stokes parameters Q (middle panel), and U (right panel). The patches are 25° on the side, centered 
at galactic coordinates (I, b) = (0°, 35°). 
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